program define estimate_IRF_LP, eclass
	
	syntax varname, [L(int 8) K(int 4) Alpha(real 0.05) ///
		min_lags_NW(int 4) Method(string) Fe_var(string) ///
		Cluster_var(string) Cluster_var_two_way(string) Time_fe_var(string)] 
	
	if "`method'" == "" local method "newey"
	
	tempvar res
	
	* Create matrices for storing results
	matrix `res' = J(`l', 7, .)
	matrix colnames `res'  = lag IRF se lb ub df_resid K
	
	* Create temporary variables containing lags and leads
	* of the variable for which the IRF is calculated
	quietly {
		forvalues ll = 1/`k' {
			tempvar `varlist'_LAG`ll'
			gen ``varlist'_LAG`ll'' = l`ll'.`varlist'
		} 
		forvalues ll = 0/`l' {
			tempvar `varlist'_LEAD`ll'
			gen ``varlist'_LEAD`ll'' = f`ll'.`varlist'
		}	
	}
	
	* Calculate local projections
	forvalues ii = 1/`l' {
	
		* Generate right-hand side variables for
		* the local projection regression
		local controls
		forvalues ll = 1/`k' {
			local controls `controls' ``varlist'_LAG`ll''
		}
		if "`time_fe_var'" != "" {
			local controls `controls' `time_fe_var'
		}
		
		local ii_1 = `ii' - 1
		
		* Estimate the regression
		if "`method'" == "newey" {
			local lags_NW = max(`min_lags_NW', `ii_1')
			quietly newey ``varlist'_LEAD`ii_1'' `controls', lag(`lags_NW')
		}
		else if "`method'" == "fe" {
			if "`cluster_var'" == "" {
				quietly reghdfe ``varlist'_LEAD`ii_1'' `controls', /// 
					absorb(`fe_var') vce(robust)	
			}
			else if "`cluster_var_two_way'" == "" {
				quietly reghdfe ``varlist'_LEAD`ii_1'' `controls',  ///
					absorb(`fe_var') vce(cluster `cluster_var')
			}
			else {
				quietly reghdfe ``varlist'_LEAD`ii_1'' `controls',  ///
					absorb(`fe_var') vce(cluster `cluster_var' `cluster_var_two_way')
			}
		}
		local se = _se[``varlist'_LAG1']
		local df_r = e(df_r)
		
		* Collect the results and calculate 
		* confidence intervals at the given significance level
		local t_crit = invt(e(df_r), `alpha' / 2)
		matrix `res' [`ii', 1] = `ii'
		matrix `res' [`ii', 2] = _b[``varlist'_LAG1']
		matrix `res' [`ii', 3] = `se'
		matrix `res' [`ii', 4] = _b[``varlist'_LAG1'] + ///
			`t_crit' * `se'
		matrix `res' [`ii', 5] = _b[``varlist'_LAG1'] - ///
			`t_crit' * `se'
		matrix `res' [`ii', 6] = `df_r'
		matrix `res' [`ii', 7] = `k'
	}
	
	ereturn matrix s `res'
	
end 
